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Abstract 

TUnfold is a tool for correcting migration and background effects in high energy physics for 
multi-dimensional distributions. It is based on a least square fit with Tikhonov regularisation and 
an optional area constraint. For determining the strength of the regularisation parameter, the L- 
curve method and scans of global correlation coefficients are implemented. The algorithm supports 
background subtraction and the propagation of statistical and systematic uncertainties, in particular 
those originating from limited knowledge of the response matrix. The program is interfaced to the 
ROOT analysis framework. 
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1 Introduction 



In high energy physics, experiments are usually performed as counting experiments, where events are 
grouped into certain regions of phase-space, also called bins. However, the kinematic properties of each 
event, such as four-momenta of particles and derived quantities, are measured only at finite precision due 
to inevitable detector effects. As a consequence, events may be found in the wrong bin. Furthermore 
there is the presence of background, such that only a fraction of the events observed in a given bin 
originates from the reaction one is interested in. 

In most cases, algorithms such as GEANT [1] are used to simulate migrations imposed by detector 
effects, whereas underlying physics processes are simulated using event generators such as PYTHIA [2]. 
After tracking the generated events through the detector simulation one is able to confront the physics 
process modelled by the event generator with the background-subtracted data. 

However, often one is interested to report results such as differential cross sections, independent of 
the detector simulation. In that case, the observed event counts have to be corrected for detector effects. 
The problem may be written as 



Vi 



Xj, 1 < i < n 



(1) 



where the m bins Xj represent the true distribution, Aij is a matrix of probabilities describing the 
migrations from bin j to any of the n bins on detector level and yi is the average expected event count at 
detector level. It is important to note here that the observed event counts yi may be different from the 
average yi due to statistical fluctuations. A schematic view is given in figure 1. The situations becomes 
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Figure 1: schematic view of migration effects and statistical fluctuations 



somewhat more complicated if there is background. In that case the yt receive an additional contribution 
from background, 

m 

Vi = ^ AijXj +bi,l<i<n (2) 

where bi is the background showing up in bin i. Both the background and the matrix of probabilities 
often suffer from systematic uncertainties which have to be considered in addition to the statistical 
uncertainties. 

One may be tempted to replace jji — > yi and x j — > Xj in equations 1 or 2 and then solve for Xj , simply 
by inverting the matrix of probabilities. However, it turns out that the statistical fluctuations of the yi 
are amplified when calculating the Xj this way. Such fluctuations are often damped by imposing certain 
smoothness conditions on the Xj . This procedure is termed "regularisation" . 

The TUnfold algorithm [3], described in this paper and interfaced to the ROOT analysis package [5], 
implements a procedure to estimate the Xj using a least square method with Tikhonov regularisation 
[4] and an optional area constraint. In order to obtain best results from the least square minimisation, 
the number of degrees of freedom, n — to, has to be larger than zero. It means that the data yi have 



f 



to be measured in finer bins than are extracted by the unfolding procedure. This condition n > to is in 
contrast to some other commonly used unfolding methods, where often the restriction n = to is imposed 
[6, 7]. Examples of unfolding algorithms which do not have the restriction n — m are [8, 9]. 

No attempt is made here to give a complete overview of the commonly used unfolding algorithms. The 
TUnfold algorithm [3] presented here compares best to algorithms based on matrix inversion or singular 
value decomposition, like [6, 10]. Alternative approaches are often based on iterative methods or on the 
use of Bayes' theorem, for example [7, 8, 9]. Many reviews on the topic can be found in literature, only 
two examples are given here [11, 12]. 



2 The TUnfold algorithm 

2.1 Definitions 

The TUnfold algorithm gives an estimator of a set of truth parameters, using a single measurement of 
a set of observables. The observables are described by a vector 1 of random variables, y. The random 
variables y are taken to have a multivariant Gaussian distribution with mean y = Ax, where x is a vector 
corresponding to the set of of truth parameters and A is a matrix. The covariance matrix of y is V yy . 
The algorithm only works if the dimension of x is less or equal to the dimension of y. Furthermore, V yy 
has to have full rank and the rows of A shall be linear independent. The algorithm returns an estimator 
x of the truth parameters x, given an observation y. The estimator x, when considered as a random 
variable, has a covariance matrix which is also calculated. It is labelled V xx . 

2.2 Algorithm 

The unfolding algorithm, as implemented in TUnfold, determines the stationary point of the "Lagrangian" 



C(x, A) =Ci +C2+C3 where (3) 

£i=(l/-Aa:) T V w - 1 (y-Aa:) ) (4) 

£2 =r 2 (x - f b x ) T (L T L)(x - f b x ), (5) 

£3 =A(y - e T x) and (6) 

Y=*Zvi, (7) 

', X (8) 

i 



The term C\ is what one expects from a least square minimisation. The vector y has n rows. The 
covariance matrix V yy of y is diagonal in many cases, such that the diagonals hold the squares of the 
uncertainties. TUnfold also supports the use of non-diagonal V yy . The vector x corresponds to the 
unfolding result and has m rows. The elements Aij of A describe for each row j of x the probabilities to 
migrate to bin i of y. The matrix A often is determined using Monte Carlo simulations. 

The term £ 2 describes the regularisation, which damps fluctuations in x. Such fluctuations originate 
from the statistical fluctuations of y, which are amplified when determining the stationary point of 
equation 3. The parameter r 2 gives the strength of the regularisation. It is considered as a constant 

throughout this paper, matrices (M) and vectors (u) are printed in bold. Matrices or vectors without indices, written 
next to each other, are multiplied. Where needed, brackets with indices are used to refer to specific elements. The notation 
M T indicates that a matrix is transposed, its rows and columns arc swapped. The inverse of M is written as M -1 . A 
vector is treated as a matrix with only one column, such that a transposed vector has only one row. The dot product of 
two vectors v± and i>2 thus is equivalent to the matrix multiplication v± T V2. Other examples are (Aa:)j = AijXj and 

(A T )ij = Aji. 
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while determining the stationary point of C. The matrix L has n columns and ur rows, corresponding to 
tir regularisation conditions. The bias vector fbXo is composed of a normalisation factor ft, and a vector 
x . In the simplest case, one has /& = 0, n# = n and L is the unity matrix. In that case, £2 simplifies 
to T 2 ||a;|| 2 , effectively suppressing large deviations of x from zero. If ft, = 1, deviations of x from Xq are 
suppressed. Choices of the matrix L different from the unity matrix are discussed in section 7. 

The term £3 is an optional area constraint. There is a Lagrangian parameter A. The sum over all 
observations is given by Y , equation 7. The efficiency vector e has m rows and is calculated from A as 
indicated in equation 8. If the area constraint is used, the normalisation of the result x, corrected for 
the efficiencies e, is thus enforced to match the total event count Y. This procedure is applied in order 
to limit possible biases on the normalisation which are present if the data y follow Poisson's statistics 
whereas the least square ansatz is strictly valid only for normal distributed measurements. The problem 
is discussed in more detail in literature, for example in [13]. 

The minimum or stationary point of C{x, A) is determined by setting the first derivatives to zero. In 
the case without area constraint, A is set to zero and only the derivatives of C\ + £2 with respect to the 
components of x are set to zero. When including the area constraint, the equations are solved for x and 
A together. The partial derivatives of £(x, A) are 

d£ g*' X) = - 2 (A T V y y _1 (?/ - Ax)) . + 2r 2 ((L T L)(x - f b x j) . - Xej, (9) 
=Y-e T x. (10) 



d£(x,X) 



dX 



(11) 



The stationary point x of £ is found as 



x = < . where (12) 

xK=o + 4Ee with area constraint 



x| A =o =E 



x\x=o without area constraint 

A 
2 

A T V y y _1 y + ^(L T L)/6^o] , (13) 
E = (A T V yy " 1 A + r 2 (L T L)) 1 and (14) 

2 = eTEe ' (15) 

In order to calculate the covariance matrix of x, given the covariance matrix of y, the corresponding 
partial derivatives are calculated 

dxk I Bf-i without area constraint 

(D y ) ki := — — =< i-(B T e ) where (16) 

oyi [Bki + (Ee)fc — e i Ee ; ' with area constraint 

B =EA T V yy ~ 1 . (17) 
The covariance matrix of the result x, originating from V yy is thus given by 

V xx = D""V yy (D x =') T . (18) 



3 Normalisation of the matrix of migrations 

In most cases, A is determined from Monte Carlo simulations. Within TUnfold, it is foreseen to initialise 
the unfolding from a matrix M of event counts, determined in a Monte Carlo event simulation, where M 
has n + 1 rows and m columns, one row more than A. The extra row is used to count those events which 
are generated in a particular bin j but are not found in any of the reconstructed bins. For the purpose 
of this paper, the extra row of M is denoted with index i = 0, whereas all other matrices and vectors 
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have indices starting from 1. In other words, the matrix elements My count the Monte Carlo events 
generated in bin j of x and reconstructed in bin i > of y, whereas the matrix elements M j count the 
Monte Carlo events generated in bin j and not reconstructed in any of the bins of y. For the unfolding 
algorithm, A and x are initialised from M as follows 

Aij —^h±^ where i > and (19) 

s j 

n 
i=0 

(jc )j =Sj. (21) 



4 Choice of the regularisation strength 

When unfolding, the strength of the regularisation, t 2 , is an unknown parameter. If r 2 is too small, the 
unfolding result often has large fluctuations and correspondingly large negative correlations of adjacent 
bins. If t 2 is too large, the result is biased towards fbX . Several methods to choose the strength of the 
regularisation are discussed in literature, for example eigenvalue analyses [14], minimisation of correlation 
coefficients [15], and the L-curve method [16]. At present, in TUnfold a simple version of the L-curve 
method is implemented to determine r 2 as well as methods to minimise global correlation coefficients. 



4.1 L-curve scan 

The idea of the L-curve method is to look at the graph of two variables L™ rve and L y nrve and locate the 
point where the curvature is maximal. These variables are defined as 

L C X UIVC =log£i and (22) 

L curve =log ^2 5 (23) 

such that L x tests the agreement of x with the data and L y tests the agreement of x with the regularisation 
condition. Note that L y urvc does not have an explicit dependence on r 2 . For t 2 — > the value of L x urve 
is minimal and L y UIVC is maximal, because £2 — > and x corresponds to the stationary point of C\ + £3. 
As t 2 gets large, L™ rvc increases whereas L™ rve is getting small, because the Lagrangian is dominated 
by £2- It is observed that the parametric plot of L y nrve against L X UIVC often shows a kink (is L-shaped). 
The kink location is chosen to determine r 2 . 

In TUnfold, the L-curve algorithm is implemented as follows: the unfolding is repeated for a number 
of points in t — log-r, thus scanning the L-curve. The curvature C of the L-curve is determined as 

d 2 L curvc dL curve _ ^curvej^curve 

C = y - 1 . (24) 

((dL™™=) 2 + (di^ u ™) 2 ) 2 

The first and second derivatives of L c x mvc (L c y mvc ) with respect to t, dL c x mvc (dL c y mvc ) and d 2 L c x urvc 
(d 2 L y UIVC ), respectively, are approximated using cubic spline parametrisations of the scan results. The 
maximum of C is finally determined with the help of a cubic spline parametrisation of C{t). 



4.2 Minimising global correlation coefficients 

A method of minimising global correlation coefficients is also implemented. Given the covariance matrix 
V xx the global correlation coefficient of a component i of x is defined as 



(y^u^h (25) 



Two sorts of correlation coefficients scans have been implemented: 
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1. minimising the average correlation: the regularisation strength t 2 is chosen such that the average 
global correlation J^. pi/n is minimised, where n is the dimension of x. 

2. minimising the maximum correlation: the regularisation strength r 2 is chosen such that the maxi- 
mum correlation maxi(pi) is minimised. 

Furthermore, it is possible to choose the covariances 

1. The covariance matrix V xx may or may not include systematic uncertainties. 

2. There is the option to partition the covariance matrix such that only parts of the matrix are used 
for the calculation of global correlation coefficients. 

3. It is possible to merge bins or groups of bins prior to calculating the pi. 

When partitioning the covariance matrix, the corresponding unused rows and columns of V xx are removed 
prior to inverting the matrix and calculating the global correlation coefficients. When merging bins of 
groups of bins, the corresponding rows or columns of the matrix are added up. 

The scan is implemented such that the unfolding is repeated for a number of points in t = logr. For 
each point the chosen correlation type (maximum or average) is calculated. The minimum is determined 
using a cubic spline interpolation. 



5 Background subtraction 

Often there is background present in the measured data y. It is worth to mention that the background 
has to include all types of events which are possibly reconstructed in one of the bins of y but do not 
originate from any of the bins of a;. In particular, part of the signal process may be generated outside 
the phase-space covered by x and thus has to be counted as background. Sometimes it is possible to 
determine background sources from the data as a part of the unfolding process, for example using a 
discriminator [17]. In order to achieve that, background normalisation factors are included as extra bins 
of the vector x, corresponding to extra columns of the matrices A, M. The background normalisation is 
then determined in the unfolding process. 

On the other hand, it is often useful to simply subtract the background prior to unfolding. Within 
TUnfold, the following method of background subtraction is implemented 



Here, the components of y° arc the data prior to background subtraction, with covariance matrix Vj y . 
The background has a normalisation factor f b with uncertainty Sf b . The background shape is described 
by a vector b and the uncertainties on the components of b are given by the vector Sb. Finally, Sij is the 
Kronecker symbol. 

The covariance matrix V yy receives contributions from the covariance matrix of y° as well as from the 
uncertainties on the background shape, the latter contributing only to the diagonal elements. In addition 
there are contributions to the covariance matrix from the background normalisation uncertainty. Because 
the background normalisation is correlated for all analysis bins, it also contributes to the off-diagonal 
elements of the matrix. 

In TUnfold, the background subtraction is generalised such that multiple background sources may be 
subtracted. The contribution of individual sources of uncertainty to the result's covariance matrix may 
be studied after unfolding. 



v =y° - f% 



(26) 
(27) 



5 



6 Systematic uncertainties on the matrix of migrations 



The matrix of migrations, A, usually receives uncertainties from various sources. First, there are statistical 
uncertainties, originating from counting the Monte Carlo events in the matrix M. Second, there may 
be systematic uncertainties, in many cases described by a variation M — > M + SM, corresponding to a 
variation of experimental conditions. 

The statistical uncertainties are bin-to-bin independent uncertainties AMy on M. They are propa- 
gated through the unfolding formalism and result in a contribution V^[' stat to the covariance matrix of 
x. Details are given in the appendix. 

A systematic variation SM is propagated to the result vector in the form of a vector of systematic shifts, 
Sx. The corresponding covariance matrix contribution is given by V^ X M = Sx(Sx) J . The calculation of 
5x is described in the appendix. TUnfold supports multiple sources of systematic variation. 



7 Choice of regular isat ion conditions 

Within TUnfold, the matrix of regularisation conditions L can be chosen with some flexibility. Three 
basic types of regularisation are supported: 

1. rows of L where only one element is non-zero, corresponding to a regularisation of the amplitude 
or size of x, 

2. rows of L where two elements are non-zero, corresponding to a regularisation of the first derivative 
of x, 

3. rows of L where three elements are non-zero, corresponding to a regularisation of the second deriva- 
tive (curvature) of x. 

The first derivatives are approximated by differences of event counts in adjacent bins, x i+ i — Xj. Similarly, 
the second derivatives are approximated by {x i+1 — x,) — (x, — Xj_i). 

When initialising TUnfold, it is possible to choose one of the three basic types of regularisation. This 
type of regularisation is then applied to all bins of x. 

1. if TUnfold is initialised to regularise on the size, L is initialised to the unity matrix. 

2. if TUnfold is initialised to regularise on the first derivatives, L has n — 1 rows and the non-zero 
elements are: Lis = — 1 and Lj^+i = 1. 

3. if TUnfold is initialised to regularise on the second derivatives, L has n — 2 rows and the non-zero 
elements are: L iti = 1, L l , i +\ = -2, L iA+2 = 1. 

On the other hand, it is also possible to choose neither of the basic types and to set up details of the 
regularisation for specific bins or groups of bins instead. 



G 



7.1 Regularisation of multi-dimensional distributions 



In many cases, x is not simply a one-dimensional distribution. Instead, the bins of x may originate from 
several distributions, for example if there are bins controlling the background normalisation in addition 
to the signal bins. Furthermore, the signal bins may originate from a multi-dimensional distribution. For 
example, the signal may have 4x3 bins in two variables Pt and 77. The vector x then has 12 bins, where 
the first 4 bins correspond to the 4 Pt bins of the first 77 bin, etc. Such a structure is not problematic 
when regularising on the size, but care has to be taken when regularising on the first or second derivatives. 

Within TUnfold there is support to initialise one-, two- or three-dimensional regularisation patterns. 
For example, when regularising the two-dimensional pattern of 4 x 3 bins from the (Pt, f]) example above 
on the second derivative, L is set up as follows: 
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(28) 



Here, rows 1-2 correspond to the regularisation of the second derivatives on Pt for the first 77 bin. 
Similarly, rows 3-4 and 5-6 act on Pt for the second and third 77 bin, respectively. Finally, rows 7-10 
correspond to the second derivatives in 77 for the four Pt bins. 



7.2 Regularisation on the density, multi-dimensional distributions 

The regularisation schemes discussed so far do not take into account the effects of non-uniform bin widths. 
Another complication arises in cases where multidimensional distributions of signal and backgrounds have 
to be mapped to the one-dimensional vectors x, y and to the matrix M. The latest version of TUnfold 
[3] addresses these issues. Multidimensional distributions are mapped on one axis of the vectors x and 
y. The regularisation conditions may be refined such that the effects of non-uniform bin widths [19] are 
taken into account. 



7.2.1 Densities 

During the unfolding, the bins of x correspond to event counts. However, often is desirable to regularise 
not on the even count but on the density. The density is calculated by dividing the number of events in 
a given bin by the width of the bin. For calculating the regularisation conditions, the number of events 
x is transformed to a density x, 

juser 

Xj — > Xj = x, x 7=r . (29) 

Ud w d 3 

The number of events Xj is divided by the multi-dimensional bin width FJ d w d j , where w^j is the bin width 
of bin j in the dimension d, as specified by the underlying multidimensional distribution. In addition, 
there is an arbitrary user function /" sor , which may be used to compensate known kinematic factors 2 . In 
TUnfold, the transformation to the density is implemented by modifying the elements of the matrix L, 



user 



where the index r is used to enumerate the ur regularisation conditions. 



L r j Y L r j x ? (^0) 



2 An example is the use of the "reduced cross section" rather than the ordinary cross section for inclusive deep-inelastic 
scattering [18]. The ordinary cross section changes by several orders of magnitude as a function of kinematic variables 
and hence is difficult to regularise. In contrast, the reduced cross section does not vary a lot, and thus is more natural to 
regularise on. 



7 



7.2.2 Derivatives 



In the case where the regularisation is made on the derivatives, the bin width may also be included in 
the approximate calculation of the derivatives. The calculation of first derivatives is modified such that 

(xj 2 — Xj ± ) — > — j (xj 2 — ajjj), (31) 

where j 2 and ji are the indices of adjacent bins of a multi-dimensional distribution and d is the dimension 
of the distribution for which the derivative is calculated. The distance between the two bin centres is Sj 2 j i 
and Ad is a normalisation constant specific to the dimension d. In TUnfold, the normalisation constant 
by default is chosen to be the average bin width in dimension d. The are relevant if derivatives 
are considered for multi-dimensional distributions, where often the derivatives along one dimension are 
different in magnitude from derivatives along another dimension. For example, in the variable Pt the 
typical bin width may be 10 [GeV], where as in r\ the typical bin width may be 0.5. In this case, the 
derivatives in P T typically are a factor of 20 smaller than those in 77, unless the normalisation A d is 
chosen appropriately. 

In analogy to the case of first derivatives, the calculation of second derivatives may be modified to 
take into account bin widths using the transformation 

(xj 3 - x j2 ) - (x j2 - x h ) -> . - — ^ ' ( 21 Zi j j (32) 

"3231 ^~ hh \ "3332 "3231 / 

where ji, an d J3 are indices corresponding to a triplet of adjacent bins of a multidimensional dis- 
tribution. The distance of bin centres and normalisation factors are defined similar to the case of first 
derivatives. 

In TUnfold, the calculation of first or second derivatives including bin widths is implemented by 
adding the appropriate modifications to the matrix L. It is possible to use the modified calculation of 
derivatives together with the density calculation explained in section 7.2.1. 



7.2.3 Example of a more complicated regularisation scheme 

Consider the use of 4 x 3 bins in (Pt, ff), where the bins borders in P t are [5, 7, 10, 15, 25] and the bin 
borders in 77 are [—2, —0.5,0.5,2]. The dimension d = 1 corresponds to P T and d — 2 corresponds to 77. 
In the example, the bin widths along p T are [2, 3, 5, 10] and those along 77 are [1.5, 1, 1.5]. The first four 
components of the vector x hold the four bins in P T of the first 77 bin, etc. The bin widths are thus given 

by 



Wl,l = W\,5 = Wi,9 = 2 
Wl,2 = wis = Wl,10 = 3 
Wl,3 = Wl,7 = Wi t u = 5 
WiA = Wl,8 = 101,12 = 10 



W2,l = W2,2 = W 2 ,3 = W 2 A = 1-5 

and w 2 , 5 = w 2 .e = ™2,7 = w 2 ,s = 1 (33) 

W2,9 = l«2,10 = 102,11 = W 2 ,12 = 1.5 



and average bin sizes are Ai = 5 and A 2 = 1.33. The distances of the bin centres are [2.5,4, 7.5] along 
Pt and [1.25, 1.25] along 77, respectively, so 



4i = 


sh = 


s 1 - 

"10,9 — 


2.5 


4,2 = 


- 51, = 


- "11.10 _ 


= 4 and 


Sl,s = 


Sir = 


"12,11 — 


7.5 



^,l=*6 1 1 2 = <J?,3 = <5f,4 = l-25 
^2 _ cl _ r2 _ r2 _ i or (^l 
"9,5 — "10,6 — "11,7 — "12,8 — 1 - zo ■ 



The resulting matrix L for the case of curvature regularisation on the density, including bin width 
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effects, then looks like 
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where the numbers have been rounded to two digits. 



8 Structure of the TUnfold software package 

TUnfold is implemented in the programming language C++ and is interfaced to the ROOT analysis 
framework. The package is organised in four classes 

TUnfold implements the basic unfolding algorithm and L-curve scan. 

TUnfoldSys inherits from the TUnfold class and adds functionality to perform background subtraction 
and propagation of systematic uncertainties. 

TUnfoldDensity inherits from the TUnfoldSys class. It adds a method to perform scans of global cor- 
relations. More important, it provides support for multidimensional binning schemes, implemented 
with the help of the class TUnfoldBinning. 

TUnfoldBinning is a class to set up binning schemes. The binning schemes are organised in tree- 
like structures. The nodes of the tree correspond to distinct channels. Each channel may hold a 
multidimensional distribution in some variables. An example of a binning scheme for the vector x 
with signal and background bins is shown in figure 2. 




Figure 2: example binning scheme with three nodes. The "generator" node is the root node. It has 
two child nodes, "signal" and "background". The "signal" node has a two-dimensional binning in two 
variables, pt and eta, whereas the background node has unconnected bins corresponding to various back- 
ground sources. 



9 Summary 

The mathematical foundations of the TUnfold software package have been presented. TUnfold can be 
used to correct measurements for migration effects using the well known mathematical techniques of least- 
square fitting and Tikhonov regularisation. For choosing the strength of the regularisation parameter, two 
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types of scanning methods are implemented: the L curve method and a flexible minimisation procedure of 
correlation coefficients. The package offers the possibility to set up non-trivial regularisation schemes for 
unfolding multi-dimensional distribution. Standard methods to subtract background and to propagate 
systematic uncertainties are also implemented. 



A Partial derivatives used for the propagation of uncertainties 



The partial derivatives of Aij with respect to Mkj are 

dM kj Sj 

The partial derivatives of x with respect to the matrix elements A^ are given by 



(36) 



dx k 




dA^ 


Ckj 


\E kj 
[E kj - 


Zi 


f( v yy~ 
\( v yy~ 



(D* y ) kl x 3 where (37) 



without area constraint 
( Ee )j( Ee )fc w ith area constraint 

e 1 Ee 



and (38) 

at 

(39) 

In order to derive this result, the partial derivatives of E with respect to the elements of the inverse E _1 
are expressed by the elements of E, 

-^- = -E ik Etj. (40) 
The partial derivative of x with respect to the regularisation parameter r 2 is 



Q Xk J ^E(L t L)(/;,:eo — x)J without area constraint 



d(r 2 ) I (E(L T L)(f b x ~ a)) ~ eTE t LT ^^'°-') (Ee) t with area constraint. 



B Propagation of systematic uncertainties 

Correlated systematic shifts are propagated in the form of systematic shifts of the result. Given a shift 
£M to the matrix M, one finds the corresponding shift dA of A using equation 19. The resulting shift on 
x is then given by 

Sx = E = C (* A ) T * - D xy (<5A)x, (42) 

Statistical uncertainties AM i3 of the elements of M may also be relevant. The calculation could be 
done by repeated application of equations 36 and 42 for each independent source of uncertainty AMy. 
However, the required computing costs are large. In TUnfold, the computation is factorised such that 
the computing cost is C(n 3 ) 

(V M, s ta t) . . = £ F . kFjkPk + Cik c jk ]T Q lk zf + 'KHirj. E Q**l 

k k i k i (43) 

- (FG T + GF T ) y - - (D xy H T + H(D xy ) T )y where 
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Q 



( ) and p 3 = Qij-> 

k k i 




(44) 



(45) 



(46) 



k OAk i 



(47) 



fe 
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